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OO ! Abstract 

Quantification of relations between measured variables of interest by statistical measures of depen- 
dence is a common step in analysis of climate data. The term "connectivity" is used in the network 

■ context including the study of complex coupled dynamical systems. The choice of dependence mea- 
sure is key for the results of the subsequent analysis and interpretation. The use of linear Pearson's 
correlation coefficient is widespread and convenient. On the other side, as the climate is widely acknowl- 
edged to be a nonlinear system, nonlinear connectivity quantification methods, such as those based on 
information-theoretical concepts, are increasingly used for this purpose. 

' In this paper we outline an approach that enables well informed choice of connectivity method for 

a given type of data, improving the subsequent interpretation of the results. The presented multi- 
step approach includes statistical testing, quantification of the specific non-linear contribution to the 
^ ■ interaction information, localization of nodes with strongest nonlinear contribution and assessment of the 

00 . role of specific temporal patterns, including signal nonstationarities. In detail we study the consequences 

00 ■ of the choice of a general nonlinear connectivity measure, namely mutual information, focusing on its 

relevance and potential alterations in the discovered dependence structure. 

We document the method by applying it on monthly mean temperature data from the NCEP/NCAR 
reanalysis dataset as well as the ERA dataset. We have been able to identify main sources of observed 
non-linearity in inter-node couplings. Detailed analysis suggested an important role of several sources 

■ of nonstationarity within the climate data. The quantitative role of genuine nonlinear coupling at this 
scale has proven to be almost negligible, providing quantitative support for the use of linear methods for 
this type of data. 
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H ! 1 Introduction 
^ . 

In climate dynamics research, analysis of time series data has a central position. Detection and quantification 
of dependence between measured or modeled variables is often of interest. 

Apart from the dependences among different physical quantities, in many contexts the task is given to 
assess the dependence among the measurements of the same physical variable (i.e. surface air temperature, 
SAT) measured at many different geographical locations. 

The motivation for such a procedure commonly stems from the need to reduce the dimensionality 
of high-dimensional o riginal data, such as in the application of Empirical Orthogonal Function analysis 



([Hannachi et all ■ l2007^ to uncover the basic modes of dynamics of climate system. On the other side, depen- 



dence quantification mi ght be used to uncov er the complex structure of the climate system using approaches 



such as graph theory (ITsonis and Roebben . 12004) . Other applications exist including those combining the 
above named, see e.g. ( Tsonis et a 
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There is a wide range of methods available for detection of dependence between variables. The most widely 
known and used is Pearson's correlation coefhcient, a measure particularly sensitive to linear dependence. 
While the Pearson's correlation detects dependence reliably in the case of multivariate Gaussian probability 
distributions, it may be suboptimal in the case of complex non- Gaussian dependence patterns. For particular 
dependence patterns (or bivariate probability distributions) , it may also fail to detect statistical dependence 
between variables of interest completely. 

Note that in the following we use the terms 'linear' and 'Gaussian' dependence interchangeably to denote 
patterns of dependence corresponding to bivariate normal distribution. While the latter term is more precise, 
the former is more commonly used in the general community together with the distinction between linear 
and nonlinear methods. 

However, alternative measures exist that are able to be tter reflect poten tial non-linear dependences. 
These include the Spearman's ordinal correlation coefhcient (|Spearmanl . [l90J) and Kendal's tau (Kendall, 



1938I ). that are designed to be sensitive to any monotonous dependence pattern, without the restriction to 



linear relationships. An ul timate altern ative to the Pearson's correlation coefficient then lies in the utilization 
of mutual information (,Shannonl . 19481 ). an information-theory based measure that is in principle sensitive 
to any dependence between variables. For this generality, mutual information is widely used to quantify 
statisti cal dependence in complex systems, and has been also introduced to the analysis of climate time 
series (|Palus and Novotnal [20o3 iDiks and Mudelseel . I2OO0I) . 

As the climatic system is highly nonlinear, it seems well m otivated to use nonl inear dependence measures 
for analysis of the measured time series, as suggested e.g. in ( Donees et al . 2009f) . This may in theory allow 
more sensitive detection and quantification of dependences, potentially uncovering new climatic phenomena. 
On the other side, nonlinear measures such as mutual information may have downsides including more 
difficult implementation and interpretation, increased computational demands and numerical stability issues. 

These considerations motivate the central question of the current report: Does the non-linear component 
of the climate time series dependences sufficiently motivate the use of nonlinear dependence measures? 

It is important to note, that the answer to this question might be complex, and certainly would be domain 
specific. To deal with this complication, we outline here first a generally applicable framework, and then 
show the results obtained by analyzing in detail a specifi c dataset of particu l ar interest. This is the monthly 
SAT data from the N CEP/NCAR r eanal ysis dataset (iKistler et aj I2O OI': 'Kal nav et all . Il99d ). as well as 



concatenated ERA-40 (jUppala et all . l2005l ) and ERA-INTERIM (|Dee et al, 2011) data. Note that this data 
has been analyzed in many recent studies, utilizing both linear and nonlinear methods, and so constitutes a 
well motivated timely and relevant example of application of this framework to guide the decision regarding 
the method choice. 



2 Materials and Methods 



Data from the NCEP/NCAR reanalysis dataset ([Kistler et all . 120011 ) have been used. In particular, we utilize 
the time series Xi{t) of the monthly mean SAT from January 1948 to December 2007 (T = 720 time points), 
sampled at latitudes Xi and longitude forming a regular grid with a step of AA = = 2.5°. The points 
located at the globe poles have been removed, giving a total of = 10224 spatial sampling points. 

To expl ore the generalizab ility of the results, analogous analysis has been c arried out for t he ERA- 
40 dataset (|Uppala et all . l2005l ) concatenated with the ERA-INTERIM dataset (|Dee et ai I2OIII ) (further 
referred to together as the ERA data) . 

To minimize the bias introduced by periodic changes in the solar input, the mean annual cycle is removed 
from the data to produce so-called anomaly time series. 

We discuss two dependence measures throughout the paper, Pearson's correlation coefficient (linear 
correlation) and (nonlinear) mutual information. 

Given two real random variables X,Y, the well-known Pearson's correlation coefficient is defined as 



Px,Y 



I E[{X - E{X)){Y - E{Y)] 
E{{X - E{X)r)E{{Y - E{Y)fy 



2 



where E{) denotes the expected value operator. The corresponding finite sample estimate is denoted by r. 
For two discrete random variables X, Y with sets of values S and T, the mutual information is defined 

as 

/(X. r, = i: log 

where ^(x) — Pr{X — a;},x £ S is the probability distribution function of X, — Pr{Y — y},y £ T 
is the probability distribution function of X and p{x,y) = Pr{{X,Y) = {x,y)},x £ S,?/ e T is the joint 
probability distribution function of X and Y. For continuous variables, mutual information is defined by 
the respective integral. However, in practice the mutual information is estimated using discretization of the 
theoretically continuous variables. 

When the discrete variables X, Y are obtained from continuous variables on a continuous probability 
space, then the mutual information depends on a partition ^ chosen to discretize the space. A c ommon 
choice is a simple box-counting algorithm based on marginal equiquantization method (jPalus et a 



i.e., a partition is generated adaptively in one dimension (for each variable) so that the marginal bins become 
equiprobable. This means that there is approximately the same number of data points in each marg inal bin. 



In th is paper we use a simple pragmatic choice of Q = 8 bins for each marginal variable (jPalus and V cimclka. 
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Mutual information is a non-negative quantity, with /(X, Y) = corresponding to independence of the 
variables X and Y, and units depending on the base of the logarithm (base 2 corresponds to 'bits', while 
natural logarithm with base e corresponds to 'nats', used here). As the estimation of mutual information 
from finite sample size is prone to sample size dependent bias, to allow quantit ative compar i son w e carry 



out an approximate correction by recalibration procedure described elsewhere (jHlinka et all . 120111 ). This 
procedure is in general based on comparison with samples of the same size and coming from populations 
with analytically established mutual information. 

To elucidate our nonlinearity assessment strategy, we first point out that for a bivariate Gaussian dis- 
tribution, the correlation p{X, Y) of the variables X, Y uniquely defines the mutual information between 
them, which is given by I{X,Y) = Ic{p{X,Y)) = — ilog(l — p^{X,Y)). However, for a general non- 
Gaussian bivariate distribution, this equation may not hold. Two cases of bivariate non-Gaussianity can be 
distinguished. 

Firstly, the 'simpler' nonlinearity consists in non-linear rescaling of one or both of the variables. Such a 
rescaling does not affect the mutual information between the variables, however, the correlation may change 
substantially. Rescaling of this type can be suspected in data e.g. due to non-linear properties of the 
measurement scale, and may be considered as a bias in the correlation estimation. A remedy commonly 
adopted is the use of Spearman rank correlation coefficient. Alternative procedure lies in preprocessing the 
data by applying a monotonous transformation to each variable separately that would render it Gaussian 
("marginal normalization"), and computing the correlation on the transformed data. 

A second, more 'substantial' type of non-Gaussianity lies in that some bivariate distributions differ from 
bivariate Gaussian not only in their marginal distributions, but also in the form of the interdependence, 
which cannot be resolved by univariate rescaling. This substantial non-Gaussianity is the key motivation 
for the use on nonlinear dependence measures, as the dependence pattern can not be recovered by only 
considering ranks or other rescaled version of the variables. 



Recently, a quantification method for such deviation from Gaussian dependence has been proposed (jHlinka et al 



201l[ ). building on the fact that for univariately Gaussian random variables X, Y, the correlation gives a lower 



bound on the mutual information by 

/(x,y) >/G(p(x,r)), 

with the minimum obtained for bivariate Gaussian distribution. In particular, one can define the neglected 
('extra- normal' or 'non-Gaussian') information 

Ie{X,Y)=I{X,Y)-Ig{X,Y)>0. 

Our investigation of the impact of nonlinear contributions to the climate dependence network is carried 
out in several steps with increasing level of detail. 
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As a first step, for each pair of local time series Xi,Xj; i,j S 1,...,A'^, the correlation r(xi,Xj) and 
mutual information I{xi, Xj) of each pair of local time series is computed, and their overall relation shown in 
a scatter plot for visual inspection of systematic deviation from the relation valid under Gaussianity. Minor 
noisy deviations are of course expected due to estimation of the quantities from finite size samples. The 
step is repeated for univariately normalized variables to isolate the effect of simple rescaling; the normalized 
variables are further used in the subsequent analysis. 

Notably, the deviation of the relation among the two estimated quantities r{xi,Xj) and I{xi,Xj) from the 
theoretical prediction under Gaussianity {I{xi,Xj) = lEir{xi,Xj))) could be still attributed to the nonlinear 
properties of the bi-variate dependences, but also to different estimator properties of correlation and mutual 
information (including residual mutual information estimator bias uncorrected by the procedure mentioned 
above) . 

To isolate the genuine nonlinearity from apparent differences due to different estimator properties, and 
to allow more robust quantitative comparison, the mutual information estimates from data X are further 
compared to mutual information estimates obtained for the respective pairs of linear surrogate data S{X). 
The surrogate data conserve the linear structure (covariance and autocovariance) and hence also the corre- 
lations of the original data, but remove non- Gaussianity (non-linearity) from the multivariate distribution. 
Thus, in the large sample size limit, the mutual information between pairs of variables in the surrogate data 
should be 

1(3,, S,) - laipiS^, S,)) = /g(p(X„X,)) < IiX,,Xj). 

Technic ally, linear surrogate data are convenien tly constructed as multivariate Fourier transform (FT) 
surrogates (jPrichard and Theileii . ll994 : |Paluslll997l) : i.e. obtained by computing the Fourier transform of the 
series, keeping unchanged the magnitudes of the Fourier coefficients (the amplitude spectrum), but adding 
the same random number to the phases of coefficients of the same frequency bin; the inverse FT into the 
time domain is then performed. 

Thus, instead of comparing the estimators of two different quantities r{xi, xj) and /(xi, xj) (or using the 
rescaled version of the first one, i.e. lEir{xi,Xj)), as an estimate of linear information based on correlation 
estimator r), we can compare the values I{xi,Xj) and I{si,Sj) obtained using the same estimator on both 
the original dataset X and its linearized (surrogate) version S. 

As an added value, generation of the surrogates allows direct statistical testing of the Gaussianity of 
the studied process, as they provide random samples with predefined covariance structure. For this second 
purpose, a set of Ngurr = 99 such surrogate datasets is generated. 

For each pair of locations i,j G 1, one can test the hypothesis that the time series Xi,Xj come from 
a bivariate linear stochastic process with Gaussian dependence among the variables by comparing the ob- 
tained mutual information estimate I{xi, Xj) with the empirical distribution of Gaussian mutual information 
estimate I{si^kT Sj^k)',k = 1 • • • , A^surr- The one-sided hypothesis is rejected at significance level a = 0.05 
if the data mutual information is higher than at least 95 (out of 99) of the surrogate mutual information 
values. 

Apart from testing statistical significance of the deviations from linearity, we also describe the strength 
and localization of the effect. This can be conveniently vizualized by estimating for each location the average 
mutual information between the location i and all other locations: for original data by J(i) = I{xi,Xj)/N 
and for a surrogate dataset by Isurrii) = X]j ^i^i' ^j)/^- Geographical rendering of the difference Isii) = 
— Isurr{i) and the relative difference lE{i)/ allows effective visual inspection of the most substantial 
nonlinearity localization. 

Such automatic localization of substantially nonlinear dependencies can be followed by inspection of the 
corresponding temporal patterns and expert assessment of their relevance for the studied phenomena. Based 
on this investigation, the use of nonlinear measures for (a subset of) the data may be recommended to exploit 
the nonlinear information in the data, or alternatively further postprocessing may be proposed to clean the 
data from spurious sources of apparent nonlinearity, as shown below. 

In particular, the intermediate results of the investigation motivated the introduction of two additional 
data preprocessing steps adopted in later sections of the paper. The first is the removal of seasonality in 
variance (variance normalization) , which removes the differences in local temperature variability in different 
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Figure 1: Relation between mutual information and correlation estimates. Each dot corresponds to one time 
series pair. For visualization purposes, a random selection of 1000 (out of over 5 million) time series pairs 
shown. 

parts of year. To achieve that, standard deviation of temperature anomalies is computed for each month 
separately and the anomaly data from given month are divided by this standard deviation. The second 
additional preprocessing step was the removal of slow trends from the data. For simplicity, only a linear 
trend was considered here. 

3 Results 

The general relation between the mutual information and linear correlation for pairs of SAT anomalies is 
visualized in Figure [TJ The mutual information is in general very strongly related to the correlation. The 
relation is even clearer after removing the estimator differences by representing the linear dependence in 
terms of mutual information in the surrogates, see Figure [2j The univariate Gaussianization does not change 
the correlation coefficients substantially (see Figure [3]) and does not affect much the relation between data 
and linear surrogates (see Figure S]). 

The deviation from purely Gaussian structure of the data has been tested separately for each pair of 
variables by comparison with a surrogate dataset distribution. More than 16% time series pairs showed 
significant at the p = 0.05 level. 
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Figure 2: Relation between total and linear mutual information estimates. Each dot corresponds to one time 
series pair. For visualization purposes, a random selection of 1000 (out of over 5 million) time series pairs 
shown. 
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Figure 3: Relation between correlation in original and marginally normalized data. Each dot corresponds 
to MI of one time series pair. For visualization purposes, a random selection of 1000 (out of over 5 million) 
time series pairs shown. 
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Figure 4: Relation between total and linear correlation estimates in marginally normalized data. Each dot 
corresponds to MI of one time scries pair. For visualization purposes, a random selection of 1000 (out of 
over 5 million) time series pairs shown. 
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Figure 5: Average contribution of nonlinear dependences. Left: average mutual information for each location. 
Middle: average nonlinear contribution to mutual information Ie- Right: average nonlinear contribution 
relative to total mutual information (Ie/I)- 
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Figure 6: Dependence patterns of six locations with high relative nonlinear dependence contribution. Top: 
total mutual information. Bottom: linear mutual information 



To investigate the observed deviations from linear dependences in more detail, we visualize the total 
and nonlinear contributions to dependence patterns in Figure [5l At first inspection one can see several well 
defined areas of relatively high nonlinear contribution to dependence patterns. These are in particular an 
extensive ring around the Anctarctica within the Southern Ocean, a few locations close to the North Pole 
(Barentz Sea, Bering Sea, Baffin Bay, Greenland Sea) and areas in Brazil and Southwest Asia. 

To understand the nonlinear dependence pattern in more detail, we select the locations with the highest 
relative nonlinear dependences and visualize both their linear and nonlinear dependence pattern with respect 
to all other locations, see Figure [S] For most of these areas, the nonlinear dependences are not generally 
stronger, but rather include additional distant locations, in contrast with the mostly local character of linear 
dependence patterns. 

This might suggest the existence of long-range i nteractions or "te leconnections" of highly or predom- 
inantly non-linear character, as discussed e.g. in ( Hsieh et al 2006h . To elucidate the nature of these 
long-range connections we inspected the bivariate distributions and time series of the variables. A represen- 
tative example is shown in Figure [71 The shape of the bivariate distribution together with close inspection 
of the time series suggests that the non-Gaussianity might be related to seasonal variability in variance of 
the signal, which further differs between the two locations. 

In this particular case, the variability at the first location is the highest in December to February, when 
it is at its lowest at the second location and vice versa in July. Thus, the information shared by these time 
series would be explainable just by the seasonal differences of dynamics and ultimately just by variation in 
local solar infiux. 
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Figure 7: Example of apparent nonlinear coupling between remote areas. The apparent nonlinearity is 
attributable to yearly cyclicity in variance, see text for details. Top: original data anomalies, middle: uni- 
variately normalized anomalies, bottom: monthly variance normalized anomalies. Left column: scatterplots 
of time series values, right column: variances of data for each month (black: 77 N, 45 E; white: 66 S, 85 W). 
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Figure 8: Average contribution of nonlinear dependences in variance-normalized data. Left: average mutual 
information for each location. Middle: average nonlinear contribution to mutual information Ie- Right: 
average nonlinear contribution relative to total mutual information (Ie/I)- 



To test this hypotheses and control for this source of bias, we re-analyze the data after normalizing the 
seasonal variance, as described in Section [2] After this additional preprocessing step, there was a marked 
decrease in detected pairs of locations with statistically significant nonlinear contribution to temperature 
dependence (8.6% at the p = 0.05 level), however, this is still more than expected by chance. The localization 
of nodes with the strongest (non)linear dependences is shown in Figure |8l 

By comparison to Figure [5] we can see that the strongest contributors to apparent nonlinear dependences 
have been mitigated by this data cleaning step. The maxima of the relative nonlinear contributions are now 
much lower and are located almost purely in the equatorial ocean regions. 

As previously, we investigate the form of a nonlinear dependence pattern related to the strongest source 
of nonlinearity in thus preprocessed data. This is shown in Figure [SI showing typical example bivariate 
distributions and time series. The source of the observed non-Gaussianity of bivariate dependence can be 
commonly tracked down to an apparent non-stationarity of the time series. In particular in the example 
shown in Figure [9] there is a strong (almost linear) trend, which might be of interest for other reasons, but 
can be considered as spurious with respect to detection of climate interactions on the time scale considered. 

This motivates additional detrending of the time series followed by yet another replication of the analysis. 

The results suggest that the trends are indeed responsible for a major part of the yet remaining apparent 
non-linearity. In particular, the amount of detected pairs of locations with statistically significant nonlinear 
contribution to temperature dependence goes further down to 6.5% at the p = 0.05 level. Similarly, the 
average nonlinear contribution to mutual information is further substantially reduced, see Figure ITOl 



3.1 ERA dataset 

The analysis was replicated for the ERA dataset. 

The general strength and distribution of nonlinear dependence within the ERA dataset is very similar 
to the NCEP/NCAR, see Figures [iTl [HI [13] and [H The fraction of pairs of nodes with significant 
nongaussianity were at the p = 0.05 level is 16.3, 7.0 and 6.1 percent in the original anomalies, variance 
normalized data and detrended variance normalized data respectively, suggesting just a slightly weaker 
contribution of the nonlinearity in the ERA dataset, and clearly not much more than expected by chance 
for linear data in the last case. 



4 Discussion 

In the previous sections we have outlined an approach for a detailed multi-step analysis of the relevance of 
nonlinearity in the dependence structure of climate time series and the results for monthly SAT reanalysis 
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Figure 9: Example of apparent nonlinear coupling between remote areas. The apparent (nonlinear) coupling 

is attributable to slow trend in the data, sec text for details. Top: monthly variance normalized anomalies, 
bottom: after removal of linear trend. Left column: scatterplots of time series values, right column: time 
series (black: N, 95 E; gray: 3 N, 30 W) 
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Figure 10: Average contribution of nonlinear dependences in detrended variance-normalized data. Left: 
average mutual information for each location. Middle: average nonlinear contribution to mutual information 
Ie = I — Isurr- Right: relative average nonlinear contribution to mutual information (Ie/I)- 
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Figure 11: Relation between total and linear mutual information estimates for the ERA dataset. Each blue 
dot corresponds to one time series pair. For visualization purposes, a random selection of 1000 (out of over 
5 million) time series pairs shown. 



13 




60 120 180 240 300 360 60 120 180 240 300 360 60 120 180 240 300 360 



Figure 12: Average contribution of nonlinear dependences for the ERA dataset. Left: average mutual 
information for each location. Middle: average nonlinear contribution to mutual information Ie- Right: 
average nonlinear contribution relative to total mutual information (Ie/I)- 
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Figure 13: Average contribution of nonlinear dependences in variance-normalized data for the ERA dataset. 
Left: average mutual information for each location. Middle: average nonlinear contribution to mutual 
information Ie- Right: average nonlinear contribution relative to total mutual information {Ie/I)- 




Figure 14: Average contribution of nonlinear dependences in detrended variance-normalized data for the 
ERA dataset. Left: average mutual information for each location. Middle: average nonlinear contribution 
to mutual information Ie = I — Isurr- Right: relative average nonlinear contribution to mutual information 
{Ie/I). 
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data. The overall picture is in general that of negligible nonlinearities in the data; the most substantial 
apparent 'nonlinearities' are attributable to nonstationarity effects. 

Therefore, the tentative suggestion with respect to choice of dependence measure for this type of data 
would be to use a linear measure (Pearson correlation coefficient), potentially after removing data nonsta- 
tionarity by preprocessing, as some of the nonstationarities may affect also Pearson correlation estimates, 
albeit differently than the nonlinear mutual information. 

An obvious question is that of generalizability of the specific findings. Explicitly, we have confirmed 
similar results in the ERA reanalysis SAT data. 

Also in daily NCEP/NCAR data we have observed relatively negligible nonlinearity, after removing 
nonstationarities in variance as well as trends (results not shown). 

While we have only showed results for the most commonly used grid with fixed angular resolution 2.5°, 
based on inspection of the spatiotemporal structure of the data we conjecture t hat an alogous results would 
be observed for other resolutions as we ll as area-corrected grids a s used e.g. by (|Donges et ai 

The present results extend those of IPalus and Novotnal (1994) who tested for possible nonlinearity in the 
dynamics of the station (Prague- Klementinum) SAT time series and found that the dependence between the 
SAT time series {x{t)} and its lagged twin {x{t + T)} was well-explained by a linear stochastic process. This 
result about a linear character of the temporal evolution of SAT time series, as well as the results of this 
study about relations between the reanalysis SAT time series from different grid-points cannot be understood 
as arguments for a linear character of atmospheric dynamics. These results rather characterize properties 
of measurement or reanalysis data at a particularly coarse level of resolution when the data reflecting a 
spatially and temporally averaged mixture of dynamical processes on a wide range of spatial and temporal 
scales are considered. For instance, a closer look on the dynamics on specific temporal scales in temperature 
and other meteorological data has led to identification of o s cillatory phenomena with nonlinear behavior, 
exhibiting phase synchronization (|Palus and Novotnal I2004L l2006l I2009L l201lt iFeliks et al I2OIOI ) . 

Also for other variables with vastly different dynamics we can expect substantial nonlinearity in bivariate 
dependencies, especially if the measurement/model has sufficient spatial and temporal resolution. Concep- 
tually, similar analysis to that presented is warranted to be carried out before the decision for use of linear 
or nonlinear dependence measure for each substantially different dataset. The work on preparation of a 
semi- automated tool for such analysis is undergoing. 

Despite the relative sparsity of the substantially nonlinear dependence patterns, even after the two 
additional preprocessing steps we have observed more than the expected proportion of significantly nonlinear 
dependences (6% instead of expected 5%). Given the number of tests carried out (several million of location- 
pairs), this small deviation likely consists a globally significant deviation from multivariate Gaussianity, 
although the intricate interdependence of the pair tests themselves makes the estimation of the exact p- value 
for a global linearity hypothesis technically very difficult. 

From practical point of view, it may be argued that the statistical determination of the above-random 
presence of apparent nonlinearity should not play a key role in method choice. Firstly, the more detailed 
quantitative analysis has already shown that even where present, the effect is relatively weak. Secondly, even 
though for some region pairs there may be a statistical indication of nonlinear dependence, it is realistic to 
suspect (given the results for the raw data, see Figures [7] and [5]) that this apparent nonlinearity may be due 
to as yet not discussed type of nonstationarity, and visual inspection of the data is required. Indeed, several 
further spurious nonlinearities have been detected due to various uncor rected proble r ns wit hin the reanalysis 
data (some of which corresponded to known problems as described in ([Kistler et alll200l[ )). 

This leads to the consideration of the strongest nonlinear pattern obs erved in the curren t data even 
after the two additional preprocessing steps on top of those carried out in (jPonges et ai I2OO9I) . As shown 
in Figure [101 ocean areas particularly in the tropical Pacific bear still a slightly elevated non-Gaussian 
contribution to dependence patterns. An example of such dependence pattern and related time series is 
shown in Figure [121 Note that the area of strongest residual non- Gaussianity roughly corresponds to regions 
implicated in ENSO dynamics; for visual comparison we plot the area and time ser ies for the El Nino 3.4 index 
(the Nino 3.4 region is bounded by 120°W-170°W and 5°S-5°N (|Trenberthl,[l997t))- The NIN03.4 SST index 



data were downloaded from from NOAA/NWS Climate Prediction Center, (http://www.cpc.ncep.noaa.gov 
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accessed November 15th, 2012). We hypothesize that t h e observed patt e rn reflects the nonlinear char acter 
of ENSO dynamics observed e.g. by lAn and JinI (|2Q04[ ) : iHannachi et all (|2003l ): iBoucharel et all (|201ll) . 



5 Conclusions 

Quantification of dependences between variables is a common task within the study of climate or other 
complex systems. The choice of dependence measures is commonly based on some theoretical assumptions 
about the underlying system. In case of climatic time series, this may lead to the choice of nonlinear 
methods due to the nonlinear nature of the underlying physical processes. However, for various reasons 
including spatiotemporal sampling or averaging, measured or modeled data might be in fact well captured 
by linear measures, with their nonlinear counterparts potentially reducing sensitivity or introducing bias. 

We have presented a multi-step approach that allows the detailed assessment of the nonlinear contribution 
to dependence patterns in a dataset, including not only statistical testing, but also quantification, localization 
and analysis of sources of this contribution. This approach can provide rationale for the decision regarding the 
choice of suitable dependence measure for given type of data, as well as direct the analyst attention to hidden 
crucial properties of the dataset. Importantly, the presented approach is quite general. It is transferable 



to ot her geoscientific datasets, as well as to other discipli nes such as neu r oscien ce, see e.g. (|Hlinka et al 



20111 ) for an earlier application of a related approach and (jHartman et all 120111) for an example focused to 
an assessment of the nonlinearity effects on graph-theoretical network characteristics. Let us n ote that the 
construction and interpr etation of graphs representing climate networks poses further challenges ( Palus et all . 



20111: iHlinka et ail2012t ). 



For monthly SAT in the NCEP/NCAR dataset and similar data, the analysis suggests that the use of 
linear dependence methods is generally sufficient, potentially after the treatment of described nonstationarity 
sources. 

This has shown that the quantitative study of amount of nonlinearity in the data provides as a side 
product valuable hints about potentially hidden specific properties of the data, that may otherwise go 
unnoticed and bias the results, were just a single method used without detailed analysis. 
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Figure 15: Example of one of the couplings with substantial residual nonlinear contribution even after removal 
of several sources of spurious nonlinearity: between two locations in the Pacific Ocean (the dependence 
pattern still substantially deviates from bivariate Gaussianity, with the estimated total and linear mutual 
information being 0.73 and 0.45 nats respectively; see bottom left for approximate location of the points). 
The apparent (nonlinear) coupHng may be attributable to nonlinearity of ENSO dynamics. Top: monthly 
SAT anomalies after seasonal variance normalization and linear trend removal (left: scatterplots of time 
series values, right: time series - black: 5S, 165W; gray: 3S, 138W), middle: monthly SAT anomalies 
without seasonal variance normalization and linear trend removal, bottom left: depiction of the respective 
locations and the SST averaging region for obtaining the El Nino 3.4 index, bottom right: El Nino 3.4 index 
time series. See text for more details. 
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